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Overview 

The  subject  of  this  talk  is  the  numerical  solution  of  the  free  RB  equations, 

M'  =  [M,n\,  M  =  vtj  +  m, 

where  M,  Q  are  skew-symmetric  matrices  and  J  is  a  diagonal  matrix  with  positive  entries. 
M  is  the  matrix  of  body  momenta 
Q  is  the  matrix  of  body  angular  velocity 
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Often  the  above  equations  are  associated  with  the  equations  that  give  the  configuration  of  the  body 
in  the  fixed  frame, 

Q'  =  QQ,  Q  e  SO (N). 

•  The  Discrete  Moser-Veselov  description  of  the  rigid  body 

•  Integrability  of  the  discrete  algorithm 

•  Backward  error  analysis  of  the  the  DMV  algorithm 

•  Higher  order  integrable  approximations 

•  Numerical  experiments  and  comparisons  with  other  methods 

•  Explicit  methods  for  the  3x3  case 


The  Moser-Veselov  discrete  version  of  the  dynamics  of  a  Rigid  Body 

The  Lagrangian  of  the  continuous  RB  equations,  is  the  kinetic  energy, 

L  =  -tr(ClT  M)  =  -tr(— Cl2J  —  ClJCl)  =  tr(f T  JCl),  (1) 

where  we  take  into  account  that  0T  =  —Cl  and  that  the  trace  is  invariant  under  cyclic  permutations. 
Following  (Marsden,  Pekarsky  &  Shkoller  1999),  discretise  Cl  =  X~xX,  where  X  e  SO (N)  is  the 
configuration  of  the  body,  using  a  finite  difference  approximation  of  the  derivative, 

n  =  x-'x  «  i xf+1(xk+1  -  xk),  xk, xk+1  e  so (x), 

which  gives 

L  «  y.rU  -  X^XM  J  -  JXj+lXt  -  XJkXMJXJMXk). 

Due  to  the  orthogonality  of  the  Xk's  and  the  cyclicity  of  the  trace,  the  first  and  the  last  term 
cancel,  and  moreover,  we  can  write 
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L  «  Jjtr (XtJXj+1). 


◄◄ 

►► 

◄ 

► 

Back 

Close 

Up  a  scaling  factor,  this  is  precisely  the  discrete  Lagrangian  of  Moser  and  Veselov. 


Consider  the  fuctional  S(X)  determined  by 

S  =  ^tr(AVA';+1) 

k 

where  X  =  {Xk}  with  Xk  G  O (TV)  and  J  is  a  diagonal  matrix.  To  obtain  the  stationary  points  of 
S,  we  consider 

5>(AVAtT+1)  -  i^tr(Afc(A'tA'tT  -/)), 

k  k 

(where  Afc  =  is  a  Lagrange  multiplier),  and  SS  =  0  becomes 

Xk+iJ  +  Xk_i  J  =  AkXk, 

from  which,  multiplying  by  Xj  on  the  right  and  taking  into  consideration  the  symmetry  of  Ak, 

Xk+1JXj  +  AVi  JXj  =  Ak  =  ATk=  XkJXj+1  +  XkJXj_,,  (2) 

hence,  the  discrete  analogue  of  the  angular  momentum  In  space, 

mk  XkJXk_^  Xk_\J Xk  , 


is  conserved. 


In  the  body  variables,  setting  u jk  =  Xj Xk_  1  e  O(N)  and  Mk  =  Xk\mkXk_  1  =  uk  J  —  Ju>k  G 
so(iV)*  (angular  momentum  w.r.t.  the  body),  (J2|  becomes 

Mk+i  =  ivkMkuJ 
Mk  =  uj  J  —  Ju>k. 


the  discrete  Euler-Arnold  equation. 


5/100 


In  the  continuous  limit:  when  tk  =  t0  +  ke,  k  =  0, 1,  2, . . ., 

•  Xk  =  X(tk ) 

•  a ik  =  XjXk_x  «  I  -  efi(tfc), 

•  Mk  «  e(Jfi  +  OJ)  =  eM(tfc), 

letting  s  — >  0,  one  obtains  the  familiar  Euler-Arnold  equations  for  the  motions  of  the  iV-dimensional 
rigid  body, 

M'  =  [M,  ft] 

M  =  Jft  +  ftJ,  Ogso  (N). 


To  solve  the  discrete  Euler-Arnold  equations  ([3]): 

•  For  k  —  0,1,2,...,  find  uk  G  SO(iV)  such  that  Mk  =  ujJ  J  —  Jujk. 

•  Update  Mk+i  =  u>kMku>J . 

By  construction,  this  algorithm 

•  preserves  exactly  momentum  and  energy  (integrable  map) 

•  is  a  second  order  approximation  to  the  continuous  rigid  body 

•  preserves  the  standard  Poisson  structure  of  T*so(N), 

{/, 9}  =  tr (M[/m, 9m\),  f,9  e  C'°°(so(iV)), 

where  f M  =  (df/dM^). 
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Note  that 


•  Also  the  IMR  is  second  order,  preserves  all  the  integrals  of  the  continuous  rigid  body. 

•  Another  much  used  method  is  a  Lie-Poisson  integrator  of  McLachlan  and  Reich  (LP2).  For 
the  3x3  RB,  it  consists  in  splitting  the  Hamiltonian 


H  =  HX  +  H2  +  H?> 


2  2 

mf  rnf, 

- 1 - 1 - - - b 

J‘2  +  J'S  J  \  +  J3 


^3 

J\  +  J2 


and  integrating  explicitly  (a  la  Strang)  the  vector  fields  of  each  split  Hamiltonian.  The  method 
is  second  order,  explicit,  preserves  the  Poisson  structure  but  does  not  preserve  H . 


Hamiltonian  error  Hamiltonian  error 
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Error  in  the  Hamiltonian  function  H  in  the 
interval  [0, 100]  and  for  h  — 


The  components  of  the  vector  mfc  for  h  — 


777-1 

0 

-m3 

7772 

m2 

= 

m  =  M  = 

m3 

0 

— 777 1 

.  m3  . 

_  -m2 

777-1 

0 

Integrability  of  the  Moser— Veselov  equa¬ 
tion 


Recall  the  Moser-Veselov  equation, 

Mk  =  uj~1  J  -  Juk ,  Mj  =  -Mfc,  ujjutk  =  I- 

in  tandem  with  the  update  Mk+ 1  =  ukMkujJ . 

•  The  Moser-Veselov  equation  Q  has  not  a  unique  solution; 

Lemma  1  (Moser  and  Veselov)  Equation  Q  is  equivalent  to  the  factorization 

(/  -  A  Mk  -  X  2J2)  =  (c ol  +  A  J)  (04.  -  A  J) 

Rewrite  the  above  factorization  as 

(/  _  A Mk  -  A2  J2)  =  ATk  (A) Ak(X)  Ak( A)  =  (t ok  -  X  J) 
and  define  the  mapping  with  image  point  Mk+1  such  that 

(/  —  XMk+i  —  X2  J2)  =  Ak(X)Aj(X). 


By  construction, 

(/  -  A Mk+1  -  X 2J2)  =  Ak{X){I  -  X Mk  -  X2J2)A~\X), 

hence  the  map  is  an  integrable  as  it  is  an  isospectral  flow. 


(4) 
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BEA  for  DMV 

Recall  the  DMV  equations  and  the  continuous  RB  equations 

Mk+ 1  =  u)kMkul,  M'  =  [M,  fi] 
Mk  —  uj-fc  J  —  J cuk ,  M  —  £IJ  H -  , 


where  04  ~  —  M7(4). 

We  wish  to  write 


Mk+\  —  —  Mk  +  /ifiWfc,  +  h2  d2  +  /^3g?3  +  h4d^  +  •  •  • , 

and  find  the  modified  vector  field 


M1  =  [M,  ft]  +  hf2(M,  ft)  +  h2f3(M,  Cl)  +  h3f4(M,  ft)  +  •  •  •  (5) 


such  that  &h(Mk)  equals  the  solution  M{tk+ 1)  at  time  tk+i  =  to  +  (k  +  l)h  of  the  modified  vector 
field  ([5]). 

To  find  we  write 

tUfc  =  exp(— hQ0  —  h2Vt\  —  h3Q2  —  h4il3  —  /i5ft4  +  ■••),  (6) 


where  ft0,fti,ft2,  •  •  • ,  are  skew-symmetric  matrices  computed  so  that 

uj  J  —  Ju>k  =  h(n(tk)J  +  Jft  (£*)). 


(7) 


we  obtain 


h(n(tk)j  +  jn(tk))  =  h(n0j  +  m0)  +  /i2(a j  +  mx  +  ±(n§ j  -  jo2)) 

+  h3(n2j  +  jn2  +  ^ [(Og^i  +  ^1^0)5  J]  t  ~b  tOq))  +  ■  ■  ■  • 


Comparing  left  and  right-hand-sides,  it  is  trivially  observed  that  the  order -h  term  disappears  if 
f20  =  Q  (to  simplify  notation,  we  omit  the  dependence  of  0  on  tk).  In  order  to  annihilate  the 
/i2-term,  we  require  that 

fti  j  +  jn  1  +  ^(n20j  -  ^0)  =  0. 

Recall  that  M  =  f 1J  +  JVL  and  hence  M'  =  Q'J  +  JO,'.  On  the  other  hand,  M'  =  [M,Q\  = 
—  (Q2J  —  JQ2).  Hence  we  can  write 

O  =  QiJ  +  jn  1  -  ^M'  =  OxJ  +  JOi  -  +  JQ!) 
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and  the  identity  is  satisfied  by  if  and  only  if 


fit  =  (8) 

In  general,  the  algorithm  to  derive  flit  for  i  =  1,  2, . . .,  is 

1.  Find  the  coefficient  of  hl+l  in  ([tJ)  and  set  it  equal  to  zero.  This  will  give  an  equation  of  the 
type  QiJ+JQi  =  CiJ+  JC/+  [A,  J],  Note  that  the  terms  C^J  +  JCi  have  an  odd  occurrence 
of  the  QjS,  while  the  terms  of  the  type  [A,  J]  have  an  even  occurrence  of  the  A,s. 

2.  Use  the  derivatives  of  M  and  Q  to  express  the  term  [A,  J]  as  CiJ  +  JCi. 


3.  Deduce  0,  =  C*  +  Q. 


n0  =  n 
n3  = 

n2  =  inw  -  in3 

n3  =  in7"  -  ^(5n2n7  +  2  nn7n  +  sirn2) 

The  functions  O* 

Once  the  0,s  are  known,  substituting  back  in  Mk+1  =  u>jMkLUk  and  using  the  well  known  identity 

OO 

exp(X)F  exp(-X)  =  expadx  Y  —  ^  —adkx(Y), 

k= 0 

where  ad x(Y)  =  [X,  Y]  and,  recursively,  ad^-(l')  =  [JX,  ad^~1(y')],  we  find  the  expressions  for  the 
functions  in  terms  of  the  Oj_!,  Oj_ 2, . . . ,  f20> 
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1  ^  / _ ]  \j  _ ^ 

di  =  ^  ^  adnfciadf2fc2  ■  ■  ■  ad nk.M,  ku  . . .  kj  e  {0, 1, . . . ,  i  -  1}.  (9) 

j= 1  ^  k1+k2-\ - f  kj=i—j 

d2  =  |([M,n,]  +  [[M,n],n]), 

d3  =  J[M,n"]  +  l[[M,n'],n]  +  i[[M,n],n']  +  i[[[M,n],n],n]  -  |[m,o3],  (10) 

•  •  •  5 


Taylor  expansion  of  the  solution  of  the  modified  equation 

Consider 

y  =  f{y )  +  hf2(y)  +  h2f3(y)  +  ■■■, 

where  /(M)  =  [M,Q]  =  [M,  J~lM\  is  the  original  vector  field  of  the  RB  equations,  where  J  is 
a  linear  operator,  defined  such  that  JVt  =  QJ  +  JQ  =  M.  Putting  y(t )  =  M(t),  we  expand  the 
solution  of  the  above  equation  in  a  Taylor  series  and  collect  corresponding  powers  of  h, 


y(t  +  h)  =  M(t)  +  hf(M)  +  h2(h(M)  +  i/7(M)f 

+  h 3  (f3[M)  +  i(/'/2(M)  +  fif(M))  +  /)(M)  +  +  ■  ■  ■ . 


where  f  is  considered  as  a  linear  operator,  /"  as  a  bilinear  operator  and  so  on  and  so  forth.  In  our 
case, 


f(z){M)  =  [z,J-1M]  +  [M,J~1z] 

=  [z,  fi]  +  [M,  J~lz\ 
f"(zi,  z2)(M)  =  2[z1,J-1z2\, 

and,  since  /  is  quadratic,  f"  and  all  the  other  higher  derivatives  equal  zero. 

At  this  point  it  is  important  to  stress  an  important  difference  between  the  expressions  for  the 
modified  vector  field  of  (Hairer,  Lubich  &  Wanner  2002)  and  ours.  While  the  vector  field  discussed 
in  (Hairer  et  al.  2002)  is  in  Mn,  hence  the  f"  is  a  symmetric  quadratic  operator,  this  is  not  the  case 
for  our  vector  field  which  is  on  matrices,  thus 


/"(/7,  f)±f"  (/,/'/)■ 


This  non-commutative  case  is  discussed  with  more  generality  in  (Munthe-Kaas  &  Krogstad  2002). 
However,  we  observe  that  all  the  terms  containing  combinations  of  /",  /'  and  /  correspond  simply 
to  higher  derivatives  of  /.  The  mixed  terms  are  treated  instead  specifically. 

After  some  algebra,  we  have 

f2  =  d2-±f'f(M ) 

=  0, 

=  ±[M,n"  -  [o,o']  —  203], 

/4  =  rf4_^M(-)-i[(///3  +  /'/) 

=  0, 

f5  =  d5-  iJIfW  -  lCf/4  +  Uf  +  Htifsf  +  /7a)) 

=  i[M,  n<to)]  -  i[M,  [o,  O'"]]  +  £[m,  O5  -  O'OO'] 

+  i[M,  [O',  O"]]  -  i[Af,  00"0]  -  i[M,  020"  +  0"02] 

+  i[M,  [03,0']]  -  ^[M,0'20  +  00'2  +  0[0,0']0], 


◄◄ 

►► 

◄ 

► 

Back 

Close 
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(11) 


time 


The  DMV  solution  of  the  RB  equations  (dotted  line),  the  exact  solution  (solid  line)  and  the  tra¬ 
jectories  corresponding  to  the  modified  vector  fields  /  +  /i2/3  (dashed  line)  and  /  +  h2h  +  h\h 
(dash-dotted  line)  in  the  interval  [0,50]  with  h  — 


Some  important  results  about  DMV 

Theorem  2  The  DMV  is  time-reversible,  hence  f2i  —  0,  i  —  1, 2, . . .. 


Theorem  3  (Moser-Veselov)  In  the  3x3  case,  the  DMV  is  a  time-reparamtetrisation  of  the  flow 
of  the  original  vector  field  of  the  rigid  body. 

Since  the  mapping  preserves  the  underlying  Poisson  structure  and  all  the  integrals  Ft  =  Ci  of  the 
system,  it  commutes  with  all  commuting  Hamiltonian  flows  generated  by  the  FjS,  M'  =  { M ,  VF*}. 
The  nonsingular  compact  level  sets  Tc  =  r\(Fj  =  cy)  consists  of  a  finite  union  of  1-dimensional 
tori  and  on  each  torus  the  DMV  mapping  is  a  shift  along  the  trajectory  depending  on  the  integral 
quantity  H2. 

Hence,  the  DMV  solves  the  modified  equation 
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M'  —  (1  +  h2r3  +  +  •  •  •  +  fi2*T2i+i  +  •  •  ')[M,  f2], 

where  h  is  the  stepsize  of  integration  and  the  r2i+1,  for  i  —  1,  2, . . .,  are  constants. 


Introduce  the  constants 


Cj,i,j  ~  J\J32  +  J[Jl  +  ^2^3> 

Cj,i  = 

Cj  =  cJtl 

A  =  ( J\  +  J2)(J\  +  ./3)(J2  +  J3) 

H2  =  (J1  +  J2)(Jl  +  J3)(J2  +  J3)H-Cj\\m\\2. 
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Theorem  4  One  has 

r3  =  ^2  ((3  det(J)tr(J)  +  ^,2)  ||m|||  +  (3Cj  +  tr(  J2))H2), 

and 

r5  =  ^I((3tr(J4)  +  27CJ,2  +  15tr(J2)CJ  +  45det(J)tr(J))Jff22 

+  (10(7^3  +  50  det(J)tr(J)Cj  +  10  det(  J)tr(  J)tr(  J2)  +  2Cjj2tr(J2)  —  28  det(  J2))\\m\\22H2 
+  (60  det(  J2)Cj  +  3CjA  +  27det(J2)tr(J2)  +  15det(J)(C'Jj2,3  +  C'J,3,2))||m||4) . 


◄◄ 

►► 

◄ 

► 

Back 

Close 

Proof.  By  direct  computation  of  /3  and  /5. 


□ 


Higher-order  integrable  methods 

For  the  original  RB  equations,  scaling  the  initial  condition  is  equivalent  to  scaling  time. 

In  our  case,  we  know  that  DMV  is  a  time-rescaling  of  the  original  RB  equation.  Therefore  we  wish 
to  rescale  the  initial  condition  to  obtain  a  better  approximation  of  the  unsealed  original  RB. 

I.C.  DMV  New  I.C.  DMV 

h(Si(tk)  j  +  jn(tk)) 

We  perform  again  the  backward  error  analysis.  We  set  now  u  =  exp(— hCl0  —  h2Qi  -\ - )  and  solve 

for  the  QjS  as  the  skew-symmetric  matrices  that  solve 

h{l  —  f3h2  +  (t2  —  f5)/i4  +  •  •  -)(f2J  +  Jf2)  =  d)T  J  —  JCj.  (12) 
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h  =  d3-^M'"  =  -f3[M:Q}+d3-^M"' 

=  -f3 [M,  n]  +  f3  =  (-f3  +  73)  [M,  O] , 

hence,  in  order  to  have  an  order-four  scheme,  we  must  set  f3  —  0  which  corresponds  to  the  choice 


◄◄ 


►  ► 


Back 


Close 


r3  =  t3. 


After  further  computations,  one  has 

f5  =  0  <-»■  f5  =  r5  -  2  r32. 

This  value  of  f5  gives  indeed  a  method  of  order  six. 

The  new  proposed  algorithms  of  order  four  and  six  are  described  below. 

The  DMV4  algorithm: 

1.  Compute  r3  and  set  M0  :=  M(t0)h/(  1  +  h2r3). 

2.  For  k  —  0, 1, . . . ,  n  —  1, 

find  the  unique  wk  as  above  such  that  Mk  =  cukJ  —  Juk 
set  Mk+ 1  =  LVkMkuk 
end 

3.  Reconstruct  Mn  :=  Mn(l  +  h2r^)/h  sa  M(tn). 

The  DMV6  algorithm: 

1.  Compute  73,  r5  and  set  f5  =  r5  —  2r|  and  M0  :=  M(to)h/(l  +  A2r3  +  A4f5). 

2.  For  k  —  0, 1, . . . ,  n  —  1, 

find  the  unique  wk  as  above  such  that  Mk  =  cuj  J  —  Jujk 
set  Mk+ 1  =  ukMkuk 
end 


3.  Reconstruct  Mn  :=  Mn(  1  +  /i2t3  +  h4r^)/h  «  M(tn). 


Some  numerical  experiments 


We  consider  with  initial  condition 


and  matrix  J  given  as 


J  = 


m0  = 


0.4165 

0.9072 

0.0588 


0.9218  0  0 

0  0.7382  0 

0  0  0.1763 


and  compare  the  DMV  explicit  scheme  with  the  Hamiltonian-splitting  method  LP2  of  (McLachlan 
1993) 
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and  the  Implicit  Midpoint  Rule  (IMR), 
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Error  at  T=100 


Step  size  h 


Error  versus  step  size  computed  at  T  —  100  for  the  methods  LP2,  IMR,  DMV,  DMV4,  DMV6. 


Error  at  T=100 
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Floating  point  operations  versus  accuracy  (T  =  100)  for  the  methods  LP2,  IMR,  DMV,  DMV4, 
DMV6.  The  roots  of  P( A)  are  recomputed  at  each  step,  use  QR  with  pivoting,  (DMV  ss  22  LP2 
per  step). 


Floating  point  operations  versus  accuracy  (T  =  100)  for  the  methods  LP2,  IMR,  DMV,  DMV4, 
DMV6  and  RATTLE6.  The  roots  of  P( A)  are  computed  once,  use  LU  instead  of  QR  (DMV  «  19 


LP2  per  step). 


Method 

h  —  — 

U  16 

h=  ± 

h  =  1.2 

h  =  2.2 

h  =  2.5 

h  =  4 

LP2 

6.7903e-03 

5.1043e-01 

1.6055e+00 

1.7002e+00 

1.9489e+00 

1.3902e+00 

IMR 

1.5494e-04 

9.9329e-03 

1.3119e-01 

3.9514e-01 

2.5276e-01 

6.1905e-01 

DMV 

1.5014e-02 

5.9899e-01 

NaN 

NaN 

NaN 

NaN 

DMV4 

1.757e-07 

7.6167e-04 

1.0785e-01 

1.6094e+00 

5.0245e-01 

5.1624e-01 

DMV6 

1.962e-10 

1.6440e-06 

7.0269e-02 

NaN 

7.0407e-01 

1.0519e-01 
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Error  for  the  various  methods  and  selected  step  sizes 


Connections  with  matrix  Ricatti  equations 

Consider  the  matrix  equation 


M  =  XJ  —  JXT. 


(13) 


Cardoso  &  Leite  (2003)  shown  that  every  solution  of  (13)  (not  necessarily  orthogonal)  is  of  the 
form 

X  =  (M/2  +  S)J~\ 

for  some  symmetric  matrix  S. 

Furthermore,  A"  is  a  orthogonal  solution  of  (13)  if  and  only  if  S'  is  a  symmetric  solution  of  the 
Riccati  equation 

S2  +  S{M/ 2)  +  (M/2)tS  -  (M2/4  +  J2)  =  0.  (14) 


Riccati  equations  are  associated  to  symplectic  matrices.  In  our  case,  the  symplectic  matrix  is 


I  Cyiripl 

If  +  J2  is  positive  definite,  it  has  been  shown  in  (Cardoso  &  Leite  2003)  that  (mJ)  has  a  unique 
solution  S  which  is  symmetric,  positive  definite,  and  such  that  the  eigenvalues  of  W  =  M/2  +  S 
have  positive  real  parts.  This  matrix  W  is  precisely  the  same  matrix  in  Moser  &  Veselov  (1991), 
from  which  one  obtains 
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(15) 


u  =  WJ~\ 


Algorithm(Cardoso  &  Leite  2003):  Compute  A",  the  unique  solution  of  (13)  in  the  special  orthog¬ 
onal  group  SO(n). 


1.  Find  a  real  Schur  form  of  H< 
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(16) 


where  Tn  and  T22  are  block  upper-triangular  matrices  such  that  the  real  parts  of  the  spectrum 
of  Tu  are  positive  and  the  real  parts  of  the  spectrum  of  T22  are  negative  definite. 


2.  Partition  Q  accordingly, 


Then,  compute 


Q  n  Q 12 

Q  21  Q  22 

S  =  Q21Q11  ■ 


3.  Compute 


X 
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Some  computational  details 

•  Compute  real  Schur  forms  by  QR  iterations  for  eigenvalues  (Golub  Hi  van  Loan  1989) 

•  Cost:  0((2N)3)  operations  (implicit  methods  for  ODEs:  0(N3)) 


N  being  the  dimension  of  M . 


The  case  N  =  3 


In  this  case, 


it  is  possible  to  find  an  explicit  spectral  decomposition  of  Hsymp i  (without  the  QR  eigenvalue 
method) 


construct  the  real  Schur  decomposition  (16)  and  hence  X  from  the  eigenstructure  of  Hf 


sympl  • 


This  yields  an  explicit  numerical  method  for  the  reduced  RB  equations. 
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The  eigenvalues  of  the  matrix  Wsympi, 

I  lym\,\ 

are  the  solutions  of  the  quadratic  eigenvalue  problem 

P( A)  =  det(A2/  -  AM  -  J2)  =  0. 

Without  loss  of  generality,  we  assume  that  J  is  diagonal,  with  entries  J1;  J2,  J3.  Then, 


(17) 
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(18) 
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=  A6  -  A4(tr( J2)  -  ||m||2)  +  A2(Cj,2  -  H2)  -  det (J2). 
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•  Reduce  to  a  cubic  equation  (compute  the  roots  explicitely) 


Schematics I  procedure 

•  Compute  eigenvalues/eigenvectors  of  Hsymp\\ 

ReA+  >  0, 

(the  eigenvectors  need  not  be  orthogonal  and  may  be  complex).  Yj ,  Y2  E  M6x3,A±  G  M3x3. 

•  Orthogonalize  the  eigenvectors  (by  Grahm-Schmidt  or  QR), 

\Y1,Y2]  =  QR, 

so  that 

HsymplQ  =  QRAR-1 

is  the  complex  Schur  form. 

•  Reduce  to  a  real  Schur  form  by  considering  real/imaginary  part  (complex  Givens  rotation). 

•  Compute  S  =  Q2iQu,  X  =  (M/2  +  S)  J~\ 


tfsympl  [Y1  Y2]  =  [  Yi  F2  ] 
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◄◄ 

►► 

◄ 

► 

Back 

Close 

•  We  don’t  need  all  the  eigenvectors,  just  Y\.  Don’t  need  R. 

•  Avoid  complex  arithmetic  altogether. 


The  numerical  DMV  algorithm 


•  m0  I— ►  hM0  =  hm0. 

•  Compute  the  eigenvalues  of  Hsymp\  =  H(Mk )  solving  for  P{ A) 

•  For  tk  =  to  +  kh,  k  =  0, 1,  2  . . . 


0  as  in  (18). 


—  Compute  the  (real)  eigenvectors  corresponding  to  A+ 

*  Compute  3  ‘quadratic’  eigenvectors  (3  matrix  factorizations,  LU/QR,  with  pivoting). 
No  need  to  compute  explicitely  L  or  Q. 

*  Compute  the  ‘dependent’  eigenvectors. 

—  Orthogonalize  the  eigenspace 

*  By  (modified)  Grahm-Schmidt  or  QR.  Only  the  Q  factor  is  needed. 

—  Compute  S  =  Q21Q11,  uk  =  (M/2  +  S)J~l 

*  Update  Mk+i  =  uj Mkujk 

•  Rescale  mN  <—  MN/h. 
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This  algorithm  produces  an  explicit  method  that  is  about  20  —  22  times  more  expensive  than  LP2, 
the  explicit  method  of  McLachlan  and  Reich. 


Concluding  remarks 

•  Explicit  algorithms  to  solve  for  the  N  —  3  free  rigid  body 

•  The  methods  are  up  to  6th  order,  completely  integrable,  possible  to  increase  to  arbitrary  order 

•  The  cost  of  the  method  is  about  10  —  22  times  more  expensive  than  the  explicit  LP2.  The 
cheaper  versions  seem  to  be  less  stable  expecially  for  large  step-size  and  long  time  computations 

•  Reconstruction  equations:  Find  the  configuration 
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The  complexive  order  is  still  2  but  the  error  is  generally  smaller.  Backward  error  analysis  reveals 
that  some  error  terms  cancel,  however  the  components  are  reparametrized  with  3  different  time 
scales. 

To  obtain  higher  order  reconstructions,  one  can  use  interpolation  of  the  vector  field  XQ  and 
its  higher  derivatives  using  the  computed  values  of  Qk. 

•  Optimal  step-size? 
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